Universal spreading of wavepackets in disordered nonlinear systems 
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In the absence of nonlinearity all eigenmodes of a chain with disorder are spatially localized 
(Anderson localization) . The width of the eigenvalue spectrum, and the average eigenvalue spacing 
inside the localization volume, set two frequency scales. An initially localized wavepacket spreads in 
the presence of nonlinearity. Nonlinearity introduces frequency shifts, which define three different 
evolution outcomes: i) localization as a transient, with subsequent subdiffusion; ii) the absence of the 
transient, and immediate subdiffusion; iii) selftrapping of a part of the packet, and subdiffusion of 
the remainder. The subdiffusive spreading is due to a finite number of packet modes being resonant. 
This number does not change on average, and depends only on the disorder strength. Spreading 
is due to corresponding weak chaos inside the packet, which slowly heats the cold exterior. The 
second moment of the packet is increasing as t". We find a = 1/3. 

PACS numbers: 05.45-a, 05.60Cd, 63.20Pw 
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The normal modes (NM) of a one-dimensional lin- 
ear system with uncorrelated random potential are spa- 
tially localized (Anderson localization). Therefore any 
wavepacket, which is initially localized, remains localized 
for all time When nonlinearities are added, NMs in- 
teract with each other 0]. Recently, experiments were 
performed on light propagation in spatially random non- 
linear optical media ^ and on Bose-Einstein condensate 
expansions in random optical potentials [4;]. 

Numerical studies of wavepacket propagation in sev- 
eral models showed that the second moment of the 
norm/energy distribution grows subdiffusively in time as 
t" [1,0] J with a ^ 1/3. Reports on partial localization 
were published as well Further numerical conduc- 

tivity studies report Ohmic behaviour at finite energy 
densities Q. 

The aim of the present work is to clarify the mech- 
anisms of wavepacket spreading and localization. We 
study two models. The Hamiltonian of the disordered 
discrete nonlinear Schrodinger equation (DNLS) 
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with complex variables V'/- The random on-site energies 
e; are chosen uniformly from the interval [—^, ^] ■ The 
equations of motion are generated by tpi = dHo / d{itlj*). 

The Hamiltonian of the quartic Klein-Gordon chain 
(KG) of coupled anharmonic oscillators with coordinates 
ui and momenta pi 
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The equations of motion arc ii/ ~ —dTix/dui^ and e; are 
chosen uniformly from the interval |] . 

We consider a wavepacket at i: = which is given by 
a single site excitation = J/Jq with e;^ — for DNLS, 
and XI — XSijg with pi — and = 1 for KG. The 
value of X controls the energy E in the latter case. 



We will use the DNLS for theoretical considerations, 
and present numerical results for both models [10]. For 
/3 = (HI) is reduced to the linear eigenvalue problem 
XAi — ciAi — (A;+i + Ai^i). The eigenvectors Ai,j are 
the NMs, and the eigenvalues Xi, are the frequencies of 
the NMs. 

The width of the spectrum {A^} is A = + 4. The 
asymptotic spatial decay of an eigenvector is given by 
A^A - e^'/« wheref(A^) < ^(0) « lOO/W"^ is the lo- 
calization length [il|. The NM participation number 
= ^/12i^ti characterizes the spatial extend - local- 
ization volume - of the NM. It is distributed around the 
mean value P^, « 3.6^ with variance ~ (1.3^)^ T5]. The 
average spacing of eigenvalues of NMs within the range of 
a localization volume is therefore A A « A/P^. The two 
scales AA < A determine the packet evolution details in 
the presence of nonlinearity. 

The equations of motion of ([1]) in normal mode space 
read 



with the overlap integral 



lAl^j^^lAiy^jA^^^l . 
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The variables (f'l' determine the complex time-dependent 
amplitudes of the NMs. 

The nonlinear frequency shift at site Iq is SX « (3. Then 
we expect three qualitatively different regimes of spread- 
ing: i) P < AA; ii) AA < /3 < A; iii) A < (3. In case 

i) the local frequency shift is less than the average spac- 
ing between excited modes, therefore no initial resonance 
overlap of them is expected, and the dynamics may - at 
least for long times - evolve as the one for /3 = 0. In case 

ii) resonance overlap may happen immediately, and the 
packet should evolve differently. For iii) the frequency 
shift exceeds the spectrum width, therefore some renor- 
malized frequencies of NMs (or sites) may be tuned out 
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FIG. 1: (color online) m2 and P versus time in log-log plots. 
Left plots: DNLS with W = 4, P = 0.1, 1, 4.5 ((b)lue, (g)reen, 
(r)ed). Right plots: KG with W — 4 and initial energy 
E = 0.05,0.4,1.5 ((b)lue, (g)reen, (r)ed). The disorder re- 
alization is kept unchanged for each of the models. Dashed 
straight lines guide the eye for exponents 1/3 (m,2) and 1/6 
(P) respectively. 



of resonance with the NM spectrum, leading to selftrap- 
ping. The above definitions are highly qualitative, since 
localized initial conditions are subject to strong fluctu- 
ations. Yet, regime iii) is also captured by a theorem 
presented in Q, which proves, that for /3 > A the single 
site excitation can not uniformly spread over the entire 
(infinite) lattice for the DNLS case. 

We order the NMs in space by increasing value of the 
center-of-norm coordinate X^, = J^i^^ti- analyze 
normalized distributions z^, > using the second moment 
m2 = ~ i')^^^ and the participation number P — 

1/^^,0^, which measures the number of the strongest 
excited sites in z^. Here P = X^i/'^-^i'- F'or DNLS we 
follow norm density distributions z„ = \4'i^\'^ / ■ 
For KG we follow normalized energy density distributions 
Zi, = hy/^^h^ with = a^/2 -|- a;^a^/2, where ai, is 
the amplitude of the vth NM and ujI = I + {K + 2)/W. 

We systematically studied the evolution of wavepackets 
for ^ and ^ ^Ts'l . The above scenario was observed very 
clearly. Examples are shown in Fig[TJ Regime iii) yields 
selftrapping (see also Figs. 1,3 in [_/1|), therefore P does 
not grow significantly, while the second moment m2 ~ 
with a « 1/3 (red curves). Thus a part of the excita- 
tion stays highly localized 0], while another part delo- 
calizes. Note that for large /3 ^ A (or similar energy for 
the KG model) almost all the excitation is selftrapped. 
Regime ii) yields subdiffusive spreading with m2 ~ 
and P ~ t"/^ (green curves). Regime i) shows Ander- 



FIG. 2: (color online) m,2 (in arbitrary units) versus time 
in log-log plots in regime ii) and different disorder strength. 
Lower set of curves: plain integration (without dephasing); 
upper set of curves: integration with dephasing of NMs, see 
text. Dashed straight lines with exponents 1/3 (no dephasing) 
and 1/2 (dephasing) guide the eye. Left panel: DNLS, W = 
4, /3 = 3(blue); W = 7, f3 = A (green); W = 10, /3 = 6 
(red). Right panel: KG, W = 10, £; = 0.25 (blue) , W = 7, 
E = 0.3 (red) , W = A, E = QA (green). The curves are 
shifted vertically in order to give maximum overlap within 
each group. 



son localization up to some time scale r which increases 
with decreasing [3. For t < t both m2 and P are not 
changing. However for i > r a detrapping takes place, 
and the packet starts to grow with characteristics as in ii) 
(blue curves). Therefore regime i) is a transient, which 
ends at some time t, and after that regime ii) takes over. 

Partial nonlinear localization in regime iii) is explained 
by selftrapping It is due to tuning frequencies of ex- 
citations out of resonance with the NM spectrum, takes 
place irrespective of the presence of disorder and is re- 
lated to the presence of exact t-periodic spatially local- 
ized states (also coined discrete breathers) for ordered 
and disordered systems (in the latter case also 
i-quasiperiodic states exist). These exact solutions act 
as trapping centers. 

Anderson localization on finite times in regime i) is ob- 
served on potentially large time scales r, and as in iii), 
regular states act as trapping centers jl^- For t > t, 
the wavepacket trajectory finally departs away from the 
vicinity of regular orbits, and deterministic chaos sets in 
inside the localization volume, with subsequent spread- 
ing. 

The subdiffusive spreading takes place in regime i) for 
t > T, in regime ii), and for a part of the wavepacket also 
in regime iii) . The exponent a does not appear to depend 
on (3. In Fig[2]we show results for m2(t) in the respective 
regime ii) for different disorder strength. Again we find 
no visible dependence of the exponent a on W. Therefore 
the subdiffusive spreading is rather universal. 

Let each NM in the packet after some spreading to have 
norm ~ n <C 1. The packet size is then 1/n » Pjy, 
and the second moment TO2 ^ 1/n^. We can think of 
two possible mechanisms of wavepacket spreading. A NM 
with index /i in a layer of width in the cold exterior. 



which borders the packet, is either heated by the packet, 
or resonantly excited by some particular NM from a layer 
with width Pi, inside the packet. Heating here implies a 
(sub)diffusive spreading of energy. Note that the numeri- 
cal results yield subdifFusion, supporting the nonballistic 
diffusive heating mechanism. 

For heating to work, the packet modes should have 
a continuous frequency part in their temporal spectrum 
(similar to a white noise), at variance to a pure point 
spectrum. Therefore at least some NMs of the packet 
should evolve chaotically in time. Ref. @ assumed that 
all NMs in the packet are chaotic, and their phases can 
be assumed to be random. With ([3]) the heating of the 
exterior mode should evolve as i4>p, « A^0^ + Pn^^'^f{t) 
where {f{t)f{t')) = 5{t — t') ensures that f{t) has a con- 
tinuous frequency spectrum. Then the exterior NM is 
increasing its norm according to ~ 0^n^t. The 

momentary diffusion rate of the packet is given by the 
inverse time it needs to heat the exterior mode up to the 
packet level: D — 1/T ^ /3^n^. The diffusion equation 
1712 ^ Dt yields 1112 ^ j3t^/'^ [Tl|- We tested the above 
conclusions by enforcing decoherence of NM phases [TtI 
and obtain 1112 ^ t^/'^ (see Figl2]). Therefore, when the 
NMs dephase completely, the exponent a = 1/2, contra- 
dicting mimeiical observations without dephasing. Thus, 
not all NMs in the packet are chaotic, and dephasing is 
at best a partial outcome. 

Chaos is a combined result of resonances and noninte- 
grability. Let us estimate the number of resonant modes 
in the packet. A NM with l^ii/p = n will excite other 
modes in first order in p-n as 



Ai/ — Ar 



(5) 



The perturbation approach breaks down, and resonances 
set in, when R^^f^ < (3n. We perform a statistical nu- 
merical analysis. For a given NM v we obtain R^^i_i,„ — 
miuyj^i^ Collecting Ru,fj.o for many v and many dis- 

order realizations, we find the probability density distri- 
bution >V(i?i/,^o) (FigEl left plot). The main result is 
that W(x ^ 0) const ^ 0. Therefore the probabil- 
ity for a mode in the packet to be resonant is propor- 
tional to (3n. On average the number of resonant modes 
in the packet is constant, proportional to /3, and their 
fraction within the packet is ~ /3n. Since packet mode 
amplitudes fluctuate in general, averaging is meant both 
over the packet, and over suitably long time windows 
(yet short compared to the momentary inverse packet 
growth rate). We conclude, that the continuous fre- 
quency part of the dynamics of a packet mode is scaled 
down by (3n, compared to the case when all NMs would 
be chaotic. Then the exterior NM is heated according 
to « X^q>^ + P^n^/^f{t). It follows 10;, |2 ~ (3^n% 



and the rate D = 1/T 
7712 ^ Dt yields 



(3 n . The diffusion equation 




FIG. 3: (color online) Left plot: Probability densities W{x) 
of NMs being resonant (see text for details) . VK = 4, 7, 10 
(from top to bottom). Right plot: Probability densities P(d) 
of peak distances between resonant NM pairs for W = 4; 
(o)range: Xc = 0.1, (b)lack: Xc = 0.01 (see text for details). 
Inset: The eigenvectors of a resonant pair of double-peaked 
states (solid and dashed lines). 



The predicted exponent is close to the numerically ob- 
served one. 

In order to clarify the nature of resonant mode pairs, 
we studied statistical properties of resonant pairs for 
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Here Xc ~ /3n is a cutoff 



m2 ~ 



a 



1/3 
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value, which decreases the more the packet spreads (in 
our simulations we reach values Xc ~ 0.01). The corre- 
sponding distributions of the overlap integral 
appear to be invariant, with an average value / w 0.05. 
At the same time, the distributions of the frequency spac- 
ing S — |Ai/ — A^ol favour smaller values the smaller Xc 
is, with an average spacing S = 0.0026 {xc — 0.1) and 
S = 0.00026 {xc = 0.01). Both NMs from a resonant pair 
with small have a multipeak, or hybrid, structure 

(inset in FigO right plot). We evaluated the probability 
density V{d) of the peak distance d, estimated with twice 
the square root of the second moment of the distribution 
zi — AliAjj^^ i (Figl3] right plot). Its average increases 

from d = 24.8 for Xc = 0.1 to d = 32.3 for Xc = 0.01. 
With the help of semiclassical tunneling rate estimates 
the level splitting of such a hybrid pair can be es- 
timated to be (5 ~ where marks the midpoint 
between the peaks. This leads to a splitting of the order 
of S ^ e~'^/^. Therefore, an increase of the packet size 
goes along with a much weaker increase of the distance 
d, in accord with our explanation from above. 

Finally we consider the process of resonant excitation 
of an exterior mode by a mode from the packet. The 
number of packet modes in a layer of the width of the 
localization volume at the edge, which are resonant with 
a cold exterior mode, will be proportional to j3n. After 
long enough spreading j3n <C 1. On average there will 
be no mode inside the packet, which could efficiently res- 
onate with an exterior mode. Therefore, resonant growth 
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can be excluded. 

The subdifFusive spreading is universal, i.e. the ex- 
ponent a is independent of /3 and W, which are only 
affecting the prefactor in Excluding selftrapping, 

any nonzero nonlinearity strength /3 will completely de- 
localize the wavepacket and destroy Anderson localiza- 
tion. The exponent a is determined solely by the degree 
of nonlinearity, which defines the type of overlap inte- 
gral to be considered in ((S]), and by the stiffness of the 
spectrum {Ajy}. We performed fittings by analyzing 20 
runs in regime ii) with different disorder realizations. For 
each realization we fitted the exponent a, and then av- 
eraged over all computational measurements. We find 
a = 0.33 ± 0.02 for DNLS, and a = 0.33 ± 0.05 for KG. 
Therefore, the predicted universal exponent a = 1/3 ex- 
plains all available data. 

Using the above approach, we estimate the growth of 
the second moment of a wavepacket which is excited on 
top of a nonzero norm density background hq . Assuming 
the wavepacket having norm density (u + hq) we find the 
rate D ~ f3'^{n + no)"* and therefore the second moment 
of the wavepacket m,2 ^ Xjv? ~ /3''(n -t- riof't. It follows 
that as long as n ^ no holds, the wavepacket spreads 
subdiffusively m2 ~ fi'^l'^t^l'^ . But once the wavepacket 
decays such that n ~ nn, normal diffusion TO2 ^ fi^n\t 
sets in, in accord with The crossover time scales to 
infinity as no is approaching zero. 

Let us generalize our results to d-dimensional lattices 
with nonlinearity order cr > 0: 

#j = Qv-i +/3iV'irv'i - X! • c^) 



Here I denotes a d-dimensional lattice vector with inte- 
ger components, and m £ D{V) defines its set of near- 
est neighbour lattice sites. We assume that (a) all NMs 
are spatially localized (which can be obtained for strong 
enough disorder W), and (b) the property W(x ^ 0) 
const ^ holds. A wavepacket with average norm n 
per excited mode has a second moment m2 I/n^/'^. It 
follows iia] 

m2 - ipH)" , a = . (8) 

The exponent a is bounded from above by amax = 1, 
which is obtained for cr — > 0. For the above studied two- 
body interaction ct = 2 we predict a{d = 2) = 1/5 and 
a{d = 3) = 1/7. 

It is a challenging task to determine the bounds of the 
domains of validity of ([5]). They will be reached when 
yV(a; ^ 0) — > 0. Is further growth prohibited then? We 
think not, because trapping an excitation in a finite 
volume must generically lead to equipartition on finite 
times due to nonintegrability. That induces a finite 
(whatever weak) chaotic component in the dynamics, 
which will heat the cold exterior. Extending the above 
resonance scenario, higher order resonances are expected 
to persist and to yield further (though slower) spreading. 
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